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Abstract 

We present numerical simulations of the 4D Edwards Anderson Ising spin glass 
with binary couplings. Our results, in the midst of strong finite size effects, suggest 
the existence of a spin glass phase transition. We present a preliminar determination 
of critical exponents. We discuss spin glass susceptibilities, cumulants of the overlap 
and energy overlap probability distributions, finite size effects, and the behavior of 
the disorder dependent and of the integrated probability distribution. 
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1 Introduction 



The Edwards- Anderson model of spin glasses [|I| turns out, not astonishingly, to be hard to 
be understood. More surprisingly (at the start) also the Sherrington Kirkpatrick mean field 
version of the model appears to be very complex. Many unprecedented features appear: 
for example the (spin glass) phase transition survives the presence of a finite magnetic field 
0] under the AT line. The Parisi Replica Symmetry Breaking (RSB) solution appears 
to describe accurately the model in the low T broken phase (and it is believed to be the 
real solution of the model). 

The relevant question is now trying to establish how many of the features of the Parisi 
solution survive when discussing finite dimensional spin glass models. For example the 
presence of a phase transition in a finite magnetic field, that is implied by the Parisi 
solution 0], is not compatible with the point of view of the droplet model ||^, where one 
expects the transition to be removed from the action of a small magnetic field. 

As usual in a complex theoretical scenario numerical simulations try to play a role (for 
a review of recent simulations of spin glass models see [^]). Here we will report numerical 
simulations that thermalize large lattices with large number of samples: we will try to show 
many signatures hinting about the presence of a phase transition, and we will be plagued 
by large finite size effects. 

Numerical simulations of spin glasses with non-zero magnetic field have now a long 
story. Maybe the first work on the subject is the one by Sourlas |]^, ^. The numerical 
simulations of [P, [1^ (see also the criticism in and the reply [1^) use the concept of 
energy overlap, that turns out to be very relevant in this study: these simulations were 
giving a first evidence, but the computers of the time allowed too small lattices for getting 
conclusive results. Simulations of [|l^ measure constant susceptibility curves and 
detect the existence of an AT line. On the contrary the work of does not detect one, 
while |]16|, |18[ are more on the yes side. Again, the work of |Jl9|] is negative on the 



existence of a transition, while claims the difficulty of establishing a firm conclusion 
(maybe a wise approach). At last recent numerical work using a dynamic approach has 
been able to claim, again, the existence of a clear transition to a spin glass phase PT| , ^ . 

Here, as we will explain in the next sections, we will find evidence that we consider 
strongly suggestive of the existence of a phase transition, and we will present very prelimi- 
nary determinations of critical exponents. We will also find that even on large lattices (on 
current standards) finite size effects are dramatically strong, and we will discuss them in 
detail. 



2 Numerical Simulations 



We will discuss here the 4:D Edwards- Anderson spin glass system with bimodal quenched 
random couplings J = ±1. Our numerical simulations have been using the parallel temper- 
ing approach |^ , that makes a real difference in the simulations of systems with quenched 
disorder. The interested reader will find for example in the last of |^ an introduction 
to optimized Monte Carlo methods (including tempering and, more in general, the multi- 
canonical approaches). 
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2.8 
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150000 


150000 


64 


66 


0.04 


1.0 


3.6 



Table 1: Parameters of the tempered Monte Carlo runs. 



We report in table (|l|) the parameters relevant for our tempered simulations: for each 
lattice size we give the number of discarded thermalization sweeps and of the sweeps 
used for measurements, the total number of samples that we have analyzed, the number of 
temperatures used in each tempered run (i.e. the number of copies of the system updated in 
parallel at different (3 values and among which the temperature values have been swapped) , 
the temperature increment, the minimum and the maximum temperature. 

The j3 values have been chosen, as customary and reasonable, in order to keep the 
acceptance factor of the tempering [3 swap of order |. In our case we have a j3 acceptance 
ratio close to 0.8 for L = 3 that goes down to a number close to 0.6 at L = 9. 

Our runs are surely well thermalized for L < 7: for each copy {3 has visited all possible 
values at least a few times, and the system has never been stuck. L = 9 is more delicate, 
and we are not sure that the data points with lower T values are fully thermalized: we have 
repeated different trial runs, with different values of the parameters, and the one we report 
here are the ones that turn out to be better thermalized. Difference among the different 
runs were in any case minor, and a possible remanence of on thermalization effects would 
affect only minor issues that we will point out in the following. Here we will only insist on 
features that are surely representing thermal equilibrium even on the L = 9 lattice. 



3 Spin Glass Susceptibilities 

In this section we will discuss the overlap and the energy overlap susceptibilities. We will 
show the signature of a spin glass like phase transition. We will discuss the location of the 
critical temperature and the determination of the critical exponents. 

We consider two real replicas of the system with spins cr, and r^, and the local energy 
operator 

^ = E , (1) 

where the sum runs over first neighbors of the site i (on a AD hypercubic lattice of linear 
size L and volume V = L^). The overlap operator is defined as 

i 

and the energy overlap operator as 
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Figure 1: % as a function of T. Empty squares for L = 9, asterisks for L = 7, crosses for 
L = 5 and horizontal bars for L = 3. 



<1E^^T.^'^' ■ (3) 

i 

This operator plays a crucial role, since it allows to distinguish a possible trivial replica 
symmetry breaking from a non-trivial breaking. In the Mean Field RSB and coincide, 
while a non-trivial behavior of q induced by the presence of interfaces would generate 
a trivial qE- Detecting a non-trivial behavior of qE is strong evidence for a RSB like 
behavior. We will consider the probability distribution of the overlap for a given sample 
of the quenched disorder, Pj{q), and the same for the energy overlap, Pj{q)- We will call 
P{q) and P^{q) the probability distribution integrated over the quenched disorder. 
We define the overlap susceptibility as 

X^^V(E{q^)-E{qf) , (4) 

where by £'( ■ ) we denote the combined operation of a thermal average and an average 
over the quenched disorder J (in the usual notation £'(■) = (■)). We define the energy 
overlap susceptibility as 

X,, = V [E{ql) - E{qEf) . (5) 

In figure @j we plot the overlap susceptibility x as a function of T for L = 3, 5, 7 and 9, 
and in figure we plot the energy overlap susceptibility xe- The statistical errors are 
computed using a jack-knife algorithm. 
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Figure 2: As in figure (0) but for Xqe- 



Botli susceptibilities sliow a divergence in tlie low T region. The comparison of figure 
(H) with the analogous plot of reference shows a difference: we do not see a cusp 
but a clear divergence (not only for L = 7 but also for L = 9 down to T = 1.0). We 
believe this is very well explained from some non complete thermalization for the larger 
lattices of reference (this possibility is proposed and discussed in reference |]20[ itself: 



a non complete thermalization has exactly the effect of smoothing the divergence, since 
correlations on large scale cannot be created): these measurements are very delicate. As 
we have already discussed we are completely safe as far as the L = 7 lattice is concerned, 
but L = 9 was clearly our limit. It is important to stress that this is the only point of 
discrepancy with reference [^: as far as other (less sensitive) quantities are concerned 
(and even for the susceptibility on the two smaller lattice sizes) we find exactly the same 
results. Reference I^Oj does not analyze the quantities based on the energy overlap. 

As we have said the divergence of the two susceptibilities is clear. We are able to follow 
the divergence on thermalized lattices down to T = 1.0. The fact that also the energy 
overlap susceptibility diverges makes a stronger case for a RSB like spin glass phase. The 
analysis of reference |^] hints that at /i = 0.4 is close or slightly lower than 1.5. We just 
note at this point that the data of figures (|l|) and (0) are fully compatible with this value. 

One can see in figures (0) and (^ that there is a temperature region where the behavior 
is asymptotic (i.e. measurements on the L = 9 lattice are compatible with the ones taken 
on smaller lattices). We use this region to fit the susceptibility as a function of the reduced 
temperature yz^r, with 
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Figure 3: \og{xg) and best fit to a power law (continuous line) and log(xgg) and best fit 
to a power law (dashed line) versus the logarithm of the reduced temperature. 



X,,{t) ^ A^J-'' . (6) 



We will use the value Tc = 1.4 (in agreement with the results of ||2T|] and with evidence 
that will be discussed later in this note). The fits in the region of T going from 2.1 to 2.6 
(that we show in figure (|^)) give 7 = 1.97 and 7^ = 1.93. Xg has a stronger L dependence 
than XqE- We plot in figure (^ the logarithm of Xg and Xgs versus log()f:), and the best fits 
to the form (H). 

The conclusion of this analysis is that the numerical data are compatible with a diver- 
gence at Tc ~ 1.4 with an exponents 7 and 7^ close to 2. 

We have also carried through a standard finite size scaling analysis, by selecting Tc 
and the critical exponents in such a way to have curves at different lattice size collapsing 
together as well as possible. As usual for not very accurate data (as is unfortunately 
frequently the case for numerical simulations of complex systems) this analysis does not 
give unambiguous results, but only hints reasonable and preferred set of values. We use 
for Xg the leading scaling form 

Xg = Lixg{LHT-Tc)) , (7) 

and the same form for Xgs- looking at Xg we find that Tc = 1.4, ^ = 0.7 and ^ = 1.3 
give a very good fit. The same values (with % = 7) also give a very good scaling behavior 
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for XqE- 

If one would take a higher value of Tc (that could be suggested by a possible interpre- 
tation of the crossover regime we will discuss in the next section, but is discouraged by 
the dynamical data of |21|), for example Tc = 1.8, we would find a higher value of of 
the order of two, and - ~ 1. Again, the finite size scaling analysis and the study of the 
asymptotic T dependence would be consistent at this effect. 



4 Skewness and Kurtosis 

In order to qualify the probability distribution of the overlap and of the energy overlap we 
will define and analyze their kurtosis (the Binder cumulant) and their skewness. For zero- 
field deterministic and disordered statistical systems the Binder cumulant g of the order 
parameter is a very good signature of the phase transition: curves of g versus T cross at 
the critical point, since the kurtosis of the probability distribution of the order parameter 
at the critical point is an universal quantity. In the infinite volume limit g goes to zero in 
the warm phase: it goes to one in a ferromagnetic phase, and to a non-trivial function in 
the broken phase of the Parisi RSB solution of the mean field spin glass theory. 

Here, since we have a non zero magnetic field, we have to consider connected expectation 
values. We define the overlap kurtosis on a lattice of linear size L as 



3 lE[iq-Eiq)r] 3 i^f) 



2 ^E{iq-Eiq)fy 2 2 



which defines x'^g^- We define the kurtosis for the energy overlap g^ in the analogous way, 
by using qE- We define the skewness of the probability distribution of the overlap as 



E{{q-E{q)y 



siT) ^ ^ = ^ , (9) 

E{{q-E{q))y xl 

(-) 

which defines Xq^ ■ In the definition of the skewness of the energy overlap probability 
distribution one substitutes qE to q. We note now that the study of the qs probability 
distribution turns out to be very important. We plot g{T) for the four different lattice 
volumes in figure (^ . In this and in next plots errors are from sample-to-sample fluctuations 
evaluated with a jack-knife analysis. 

The difference with the usual zero field picture is strong. There is a clear change of 
regime close to T = 2 (the critical point at h = 0). The L = 3 system is small and 
different, and never really gaussian in our T range. For L > 5 g becomes non-trivial when 
T becomes smaller than 2. The g-kurtosis in this region does not change much with size. 
In the statistical error the values of the L = 5, 7 and 9 lattice are compatible (but maybe 
at T ^ 1 where also small non-equilibrium effects have to be accounted for). The fact that 
the kurtosis does not depend on L and is non-trivial is very clear from our data. Again, 
the traditional crossing behavior is completely absent (likely even in the infinite volume 
limit) in our data. Two things must be stressed here: first of all that the existence or non- 
existence of a crossing also depends from the shape of the critical asymptotic probability 
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Figure 4: As in figure (|I]) but for tlie overlap kurtosis g. 

distribution, and second tliat we cannot exclude, and on tlie contrary we clearly detect 
(see for example next figure, (^) the existence of very strong finite size effects. 

Figure (||), where we plot the energy overlap kurtosis, is dramatically and delightfully 
different from figure (^. It is already very interesting to look at the high T region: here 
only the L = 9 lattice starts to be gaussian, while smaller lattices have a strongly non- 
gaussian behavior. In the cold region again there is a strong finite size dependence (we will 
see in the next section that there is a finite size double peak structure that is disappearing 
on large lattices). Here there is a crossing, even if it is inverted as compared to usual, h = 
systems, where in the warm phase the kurtosis curves become lower with increasing lattice 
size: in our case for small sizes and high T we have a negative kurtosis, that tends to zero 
from below when the size increases. This real, inverted crossing is at T higher than 1.5, 
but we believe it should be taken cum grano salis as far as the exact determination of the 
critical point is concerned. The finite size effects make themselves clear in the non-gaussian 
behavior in the warm phase, and in a peak in the cold phase, that shrinks and shifts with 
increasing lattice size (see figure (^). A value of Tc ~ 1.5 is very compatible with this 
picture. 

The behavior of the skewness for the overlap (^) and for the energy overlap (|^) repeats 
a similar pattern. We find a symmetric behavior of the overlap in the warm phase, and 
curves collapse to a non-trivial shape in the low T regime. The energy overlap skewness, 
on the contrary, starts to become zero in the warm region only on our larger lattice, L = 9, 
and heavily depends on L in the low T region. Again a peak close to T ~ 2 is shrinking 
with increasing lattice size. 
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Figure 7: As in figure (|l]) but for tlie energy overlap skewness s^. 

In tlie next section we will discuss the full probability distribution. The analysis we 
have discussed here strongly suggests the existence of a RSB phase. We have been able to 
get hints for the value of Tc and of the critical exponents, but such estimates, that surely 
improve on existing ones, have to be considered only as hints. 



5 Probability Distributions 

Before starting a detailed discussion of the behavior of the probability distribution of the 
order parameter for a given sample Pj{q) and of the disorder averaged P{q) we briefly 
discuss finite size effects on, for example, {q)L- We have analyzed the size dependence of 
{q)L and {qE)L- We have tried a fit of the type (see Jl^ for a detailed discussion of the 
issue) 



{q)L - {q)oo + , 

{qEjL ^ {qE)oo + jj^ ■ (10) 

In the Parisi solution of the mean field theory dq = dqj^: corrections that scale as V^'^ and 
the fact that the upper critical dimension is 6 imply that in the mean field limit we expect 
dq = 2. 

We find that for the overlap the best fit to dq (using L = 3, 5, 7 and 9) increases with T 
from 2.3 ±0.5 at T = 1.0 to 3.2 ±0.2 at T = 1.6, and to a number larger than 4 for T among 
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Figure 8: P{q) at L = 9 for different T values: from the left are the curves for T = 2.4, 
2.2, 2.0, 1.6, 1.4, 1.2, 1.0. 



2 and 3 (for T > Tc we expect an exponential decay, but on a finite lattice with a finite 
number of points we can fit an effective exponent), is larger: here finite size effects are 
smaller, and the exponent more difficult to determine. In the region of T = 1.6, 1.8, 2.0 we 
find an exponent close to 4, that increases in the high T region. For example at T = 1.6, 
that asymptotically is probably marginally off-critical but very close to the estimated Tc, 
where we have a clean determination of both exponents, we have a ratio ~ | (to be 
compared with the 1 that one finds in the mean field theory). The behavior of {q)L does 
not look compatible with a pure exponential, while {qE)L would also be compatible with 
an exponential decay. We stress again that since the expected finite size corrections have 
a very complex pattern (different terms could be leading for different L values) the exact 
theoretical significance of this numerical result is not clear (see [0 for a discussion), but 
for the fact that we find that the two exponents are not very different from the ones that 
are found in mean field. For a comparison of the equilibrium value E{q), and the dynamical 
value of the overlap, qn, see figure two of ref. (our data for E{q) are very similar to 
the ones of 13). 

Let us examine in some more detail now the full P{q) and P^{q) and the probability 
distributions for a given sample, Pj{q) and Pj{q). In figure (||) we plot P(g) at L = 9 
for different T values. At T = 2.4 P is gaussian and symmetric. For lower T values we 
get a strong deviation from a Gaussian behavior. We will discuss later the fact that finite 
volume effects are very strong (and they turn out to be the most serious limitation in the 
case of the simulations we are discussing here). We can already notice that at T = 1 from 
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Figure 9: As in figure @ but for P^{q). 



the results of |^ we expect that the minimum value allowed for the overlap is Qmin — 0.55. 
On the contrary even on our larger size, L = 9, we get a long tail that basically goes down 
to g = (we will see later that it goes down to —1 for the smaller lattice sizes). 

A very similar pattern holds for P^{q) in figure (|^). Here even at L = 9 and at T = 2.4 
(deep in the warm phase) the probability distribution is not yet fully gaussian (as one can 
also see from the kurtosis and the skewness). As for P{q) the distribution becomes very 
asymmetric at low T, with a large tail towards small overlaps. 

It is interesting to analyze in more detail the size dependence of the probability distri- 
butions. In figure ([T0|) we show P{q) versus g at T = 2.0 for L = 3, 5, 7 and 9. This is the 
point of the transition in the h = model, asymptotically in the warm region for h = 0.4. 
Clearly for L = 3 the distribution is far from Gaussian: only at L = 7 one sees a Gaussian 
behavior. 

The same function at low T is very different. We show P{q) versus q at T = 1.0 in 
figure ([TTD . Also here finite size effects are very large, but P{q) does not show any sign of 
convergence to a gaussian behavior. In the mean field RSB solution one would expect a 
delta function at qmin (that, as we already said, at T = 1.0 is close to 0.55 
least for L < 9, we do not see a S function like contribution at g < qmax- Only at L 



here, at 



9 we 



see a small bump at g ~ 0.2: we see it in all our different runs, but we cannot exclude (and 
on the contrary we believe it is possible) that it is due to lack of complete thermalization 
of the L = 9 data at the lower T values (as we have discussed even if quantities like (g^) 
are well behaved and apparently thermalized we cannot completely exclude a very small 
effect of this kind). Also we will discuss a similar effect, that turns out to be due to the 
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Figure 10: P{q) versus q at T = 2.0 for L = 3, 5, 7 and 9 (curves from bottom to top). 




-1 -0.5 0.5 

q 



Figure 11: As in figure (p!0|) but T = 1.0. 
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Figure 12: As in figure (|TOD but Pe- 



finite size of the lattice, for the energy overlap. 

In figure (|^) we show Pe{<1e) versus g at T = 2.0 for L = 3, 5, 7 and 9. Again, on 
small lattices even at warm T Pe is not symmetric (the tail is in this case for large values 
of the overlap). 

In figure (|T3|) we show that in the cold region (T = 1) the energy overlap has even 
stronger finite size effects than the spin-spin overlap. Here a spurious peak at low values of 
qE is very strong for small L values (it carries the thirty percent of the weight at L = 3), 
and becomes smaller and smaller on larger lattices. 

We also show, in figures (14) and ([T5|) two individual P{q) for two different realizations 
of the quenched disorder at L = 9, T = 1. It is clear that different samples can have a very 
different equilibrium probability distribution, and the system does not look self-averaging. 
In the two examples we show, for example, it is clear that, like in mean field theory, we can 
find systems with one maximum of P{q) and systems with a complex structure of P{q), 
with many local maxima. The integrated P{q) is not originated from a trivial sum of very 
similar individual Pj{q), but from the sum of individual very different distributions. It 
is clear that this feature, crucial in the RSB mean field picture, is shared by the finite 
dimensional system in a non-zero magnetic field. 



6 Conclusions 



We have discussed accurate numerical simulations of the 4D Edwards Anderson spin glass in 
magnetic field. Our results hint very strongly the existence of a spin glass phase transition. 
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Figure 13: As in figure (|TOD but Pe and T = 1.0. 
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Figure 14: A first Pj{q) for one disorder sample, L = 9 and T = 1. 



15 



0.2 0.4 0.6 0.8 1 

q 



Figure 15: A second Pj{q) for a different disorder sample, L = 9 and T = 1. 

Susceptibilities grow strongly at low T, and fits to a divergent behavior are very good. 
Cumulants of the overlap and energy overlap probability distribution like the kurtosis and 
the skewness show a clear change of regime in the region of temperatures lower than the 
h = critical point. Finite size effects are dominated by power laws that are similar to the 
ones of the mean field theory. Probability distributions are non-trivial in the low T region, 
and, what is most important, different samples behave clearly in very different ways. The 
energy overlap follows the usual overlap in this RSB like behavior, making the possibility 
of a non-trivial behavior caused by interfaces quite unplausible. 

There are also differences with the usual RSB like picture at /i = 0. For example here 
the Binder cumulants do not cross (and the pictures of the P{q) show why). What is more 
impressive and relevant is the presence of very strong finite size effects (even on lattices of 
size = 9^, that had never before been thermalized in a numerical simulation). These 
effects are far larger than in the h = case. This is the most important limitations of the 
present numerical simulation (and of the physics conclusions one can draw from them): 
when finite size effects are as strong as we have shown is very difficult to be sure that any 
asymptotic behavior has been observed. Because of that all the quantitative results we 
quote for critical exponents and temperatures have to be taken as simple indications. The 
other real problem, as far as the coincidence with the RSB mean field picture is concerned, 
is that even for L = 9 we do not see any trace of a (5-function at q = qmin'- even if on 
theoretical grounds we expect this peak to be smaller than the one at qmax ^]^, and if we 
know that finite size effects are very strong, and if we have a small bump in P{q) at L = 9 
at low q (that we cannot take too seriously), this a worrying point, that is there to demand 
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further clarification. 
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